#ifndef SOLVERBASE_HPP
#define SOLVERBASE_HPP

#include "crtpBase.hpp"
#include "meshTypes_fwd.hpp"
#include "simulationTime.hpp"
#include "solverTiming.hpp"
#include "field.hpp"
#include "constantField.hpp"
#include "sourceField.hpp"
#include "edgeTraversal.hpp"
#include "invalidByDefault.hpp"


#include <array>
#include <experimental/memory>
#include <iostream>
#include <vector>

namespace hms
{

template<typename SolverType>
class SolverTraits {};

template<typename SolverType>
struct SolverSettings : public InvalidByDefault {
	using TraitsType = SolverTraits<SolverType>;
	using SourceFieldName = typename TraitsType::SourceFieldName;
	using FluxScheme      = typename TraitsType::FluxScheme;
	std::array<bool, static_cast<int>(SourceFieldName::NNames)>
		useSourceFields;

	FluxScheme fluxScheme;

	EdgeTraversal edgeTraversal;
};

template<typename _SolverType>
class SolverBase : public CrtpBase<_SolverType>
{
public:
	using SolverType        = _SolverType;
	using TraitsType        = SolverTraits<SolverType>;
	using SettingsType      = typename TraitsType::SettingsType;
	using MeshType          = typename TraitsType::MeshType;
	using ConstantName      = typename TraitsType::ConstantName;
	using FieldName         = typename TraitsType::FieldName;
	using ConstantFieldName = typename TraitsType::ConstantFieldName;
	using SourceFieldName   = typename TraitsType::SourceFieldName;
	using FluxScheme        = typename SettingsType::FluxScheme;

	static_assert( std::is_base_of_v<MeshBase<MeshType>, MeshType>,
		"MeshType must be derived from MeshBase!"
	);

	static constexpr int nConstants =
		static_cast<int>(ConstantName::NNames);

	static constexpr int nFields =
		static_cast<int>(FieldName::NNames);

	static constexpr int nConstantFields =
		static_cast<int>(ConstantFieldName::NNames);

	static constexpr int nSourceFields =
		static_cast<int>(SourceFieldName::NNames);

protected:
	std::array<scalar, nConstants> m_constants;
	std::array<Field<MeshType>, nFields> m_fields;
	std::array<ConstantFieldVar<MeshType>, nConstantFields> m_constantFields;
	std::array<SourceFieldVar<MeshType>, nSourceFields> m_sourceFields;
	/* even though only the domain entries of the storage field are required,
	 * giving it the same size as the mesh allows for unified indexing */
	Field<MeshType> m_storage;
	#ifdef SOLVER_STORE_FLUXES
	Field<MeshType> m_fluxes;
	#endif

	std::experimental::observer_ptr<const Time> m_time {nullptr};

	bool m_isInitialised {false};

	SettingsType m_settings;

	ExecTimings m_execTimings { ExecTimings::Zero() };

public:
	SolverBase();

	SolverBase(
		std::array<scalar, nConstants>&& constants,
		std::array<Field<MeshType>, nFields>&& fields = {},
		std::array<ConstantFieldVar<MeshType>, nConstantFields>&& constantFields = {},
		std::array<SourceFieldVar<MeshType>, nSourceFields>&& sourceFields = {}
	);

	void assertInitialisedMembers() const;

	void assertSameMesh() const;

	void checkFieldSizes() const;

	void initialise();

	bool isInitialised() const;

	scalar calculateTimeStep() const;

	void calculateFluxes();

	void calculateSources();

	void calculateNewState();

	void updateBoundaries();

	void advance( Time& time );

	auto time() const -> const Time&;

	void time( const Time& t );

	auto mesh() const -> const MeshType&;
	

	auto constants() const -> const std::array<scalar, nConstants>&;
	auto constants() -> std::array<scalar, nConstants>&;
	auto constant( ConstantName name ) -> scalar&;
	auto constant( ConstantName name ) const -> scalar;

	auto fields() const -> const std::array<Field<MeshType>, nFields>&;
	auto fields() -> std::array<Field<MeshType>, nFields>&;
	auto field( FieldName name ) const -> const Field<MeshType>&;
	auto field( FieldName name ) -> Field<MeshType>&;

	constexpr static auto zFieldName() -> FieldName;
	auto zField() const -> const Field<MeshType>&;
	auto zField() -> Field<MeshType>&;

	auto constantFields() const ->
		const std::array<ConstantFieldVar<MeshType>, nConstantFields>&;
	auto constantFields() ->
		std::array<ConstantFieldVar<MeshType>, nConstantFields>&;
	auto constantField( ConstantFieldName name ) -> ConstantFieldVar<MeshType>&;
	auto constantField( ConstantFieldName name ) const
		-> const ConstantFieldVar<MeshType>&;

	auto sourceFields() const ->
		const std::array<SourceFieldVar<MeshType>, nSourceFields>&;
	auto sourceFields() ->
		std::array<SourceFieldVar<MeshType>, nSourceFields>&;
	auto sourceField( SourceFieldName name ) -> SourceFieldVar<MeshType>&;
	auto sourceField( SourceFieldName name ) const
		-> const SourceFieldVar<MeshType>&;

	auto useSourceFields() const -> const std::array<bool, nSourceFields>&;
	void useSourceFields( std::array<bool, nSourceFields>&& );

	constexpr bool useSourceField(SourceFieldName) const;
	void useSourceField(SourceFieldName, bool);

	template<typename Visitor>
	void visitUsedSourceFields( Visitor&& ) const;
	
	template<typename Visitor>
	void visitUsedSourceFields( Visitor&& );

	void updateSourceFields();

	auto storage() const -> const Field<MeshType>&;
	auto storage() -> Field<MeshType>&;

	#ifdef SOLVER_STORE_FLUXES
	auto fluxes() const -> const Field<MeshType>&;
	auto fluxes() -> Field<MeshType>&;
	#endif

	auto settings() const -> const SettingsType&;
	auto settings() -> SettingsType&;

	auto execTimings() const -> const ExecTimings&;

	auto logger() const -> log::logger&;

private:
	void validate() const;
};


} // end namespace hms

#include "solverBase.tpp"

#endif